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ABSTRACT 

We use numerical simulations in a ACDM cosmology to model density profiles in 
a set of sixteen dark matter haloes with resolutions of up to seven million particles 
within the virial radius. These simulations allow us to follow robustly the formation 
and evolution of the central cusp over a large mass range of 10 11 to 10 14 M , down to 
approximately 0.5% of the virial radius, and from redshift 5 to the present, covering 
a larger range in parameter space than previous works. We confirm that the cusp of 
the density profile is set at redshifts of two or greater and remains remarkably stable 
to the present time, when considered in non-comoving coordinates. 

Motivated by the diversity and evolution of halo profile shapes, we fit our haloes 
to the two parameter profile, p oc ?- r ^ r , pti^ c r / r . )13--y ; where the steepness of the 
cusp is given by the asymptotic inner slope parameter, 7, and its radial extent is de- 
scribed by the concentration parameter, c 7 (with c 7 defined as the virial radius divided 
by the concentration radius). In our simulations, we find 7 ~ 1.4 — 0.08Logio(Af/M*) 
for haloes of 0.01M* to 1000M*, with a large scatter of A7 ~ ±0.3, where A/* is the 
redshift dependent characteristic mass of collapsing haloes; and c 7 ~ 8.(M/A/*)~°' 15 , 
with a large Af/AT* dependent scatter roughly equal to ±c 7 . Our redshift zero haloes 
have inner slope parameters ranging approximately from r _1 (i.e. Navarro, Frenk, & 
White) to r -1 - 5 (i.e. Moore et al. ), with a median of roughly r~ 13 . This two pa- 
rameter profile fit works well for all types haloes in our simulations, whether or not 
they show evidence of a steep asymptotic cusp. We also model a cluster in power law 
cosmologies of P oc k n , with n = (0, -1, -2, -2.7). Here we find that the concentration 
radius and the inner cusp slope are a both function of n, with larger concentration 
radii and shallower cusps for steeper power spectra. 

We have completed a thorough resolution study and find that the minimum re- 
solved radius is well described by the mean interparticle separation over a range of 
masses and redshifts. The trend of steeper and more concentrated cusps for smaller 
Af/M* haloes clearly shows that dwarf sized ACDM haloes have, on average, signif- 
icantly steeper density profiles within the inner few percent of the virial radius than 
inferred from recent observations. 

Code to reproduce this profile can be downloaded from 
[http://www.icc.dur.ac.uk ~~ reed /profile . html 
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1 INTRODUCTION 

The mass distribution of dark matter haloes provides a di- 
rect probe of the nature of the dark matter particle, as the 
inner structure of dark matter haloes is particularly sensi- 
tive to the dark matter properties. For example, warm dark 
matter should produce lower density halo cores than cold 
dark matter (CDM) because of the phase density ceiling 
introduced by the non-zero thermal velocity of warm par- 
ticles (e.g. Tremaine & Gunn 1979). The variation of peak 
halo phase density with halo mass is also dependent on the 
"coldness" of the dark matter particle (e.g. Lake 1989 and 
references therein). Spectroscopic observations of stellar mo- 
tions in galaxies, lensing properties and X-ray temperature 
maps of cluster cores each can provide a measurement of the 
central dark matter distribution in haloes, albeit with some 
uncertainty inherent in inferring the dark matter distribu- 
tion from properties of baryons or baryon dominated re- 
gions. CDM haloes in N-body simulations consistently have 
a steep central cusp where the density rises as r" 1 (Navarro, 
Frenk, & White 1996, 1997, NFW hereafter; Huss, Jain, & 
Steinmetz 1999; Power et al. 2003), r" 15 (Moore et al. 1998, 
1999, M99 hereafter; Taylor & Navarro 2001; Governato, 
Ghigna, & Moore 2001; Fukushige & Makino 2001, 2003), 
or somewhere in between (e.g. Klypin et al. 2001; Fukushige, 
Kawai, & Makino 2003; Hayashi et al. 2003; Navarro et 
al. 2004). These numerical findings appear in conflict with 
the most direct observational results. Rotation curves of low 
surface brightness (LSB) dwarfs consistently yield density 
profiles with nearly constant density cores (e.g. Flores & 
Primack 1994; Moore 1994; Salucci & Burkert 2000; de Blok 
et al. 2001). Studies of more luminous galaxies imply sim- 
ilar problems (e.g. Salucci & Burkert 2000; Salucci 2003; 
Gentile et al. 2004). The disagreement with rotation curves 
may indicate an insurmountable problem with CDM mod- 
els, may be due to uncertainties in measuring accurate stel- 
lar curves at just ~ 1% of the virial radii (van den Bosch 
et al. 2000; van den Bosch & Swaters 2001; see however 
e.g. Simon et al. 2003), or may perhaps reflect some system- 
atic bias common to all high resolution N-body simulations. 
Alternatively, the disagreement may be due to a problem 
with common assumptions made when reconstructing mass 
profiles from circular velocity data (Hayashi et al. 2003). 
Strong gravitational lensing in clusters can potentially pro- 
vide a direct measurement of the halo mass profile, and in- 
deed central mass profiles for several lensing clusters have 
been calculated (Tyson, Kochanski, & Dell'Antonio 1998; 
Shapiro & Iliev 2000; Sand, Treu, & Ellis 2002; Gavazzi 
et al. 2003; Sand et al. 2004), but have yielded conflict- 
ing results. Cluster density profiles inferred from Chandra 
luminosity-temperature mapping have steep cusps that are 
inconsistent with the flat cores observed from LSB rotation 
curves (Lewis, Buote & Stocke 2003), though this method 
is sensitive to models of the intracluster gas. In sum, many 
observational studies suggest a flatter profile than predicted 
by CDM models, but observations have not yet converged 
upon a basic shape of the density profile (e.g., Jimenez et 
al. 2003). 

NFW found that CDM haloes have a "universal" den- 
sity profile that is independent of mass, cosmological pa- 
rameters, and the initial density fluctuation spectrum with 
significant scatter from halo to halo, 



' (r/r.)(l + r/r.)a' w 

where r s and p s is a characteristic inner radius (the concen- 
tration radius) and inner density, respectively. Fukushige & 
Makino (1997), based on a single halo with ~ 10 6 particles, 
as opposed to the ~ 10 4 particles of the NFW study, found 
a profile with slope between r _1 and r~ 2 . M99, using results 
from a series of six ~ 10 6 particle haloes, also found a profile 
steeper than r , and proposed the following profile: 



' (r/r s )^[l + (r/r s )^Y W 

The NFW and M99 profiles are both specific cases of a three 
parameter profile family proposed by Hernquist (1990), and 
further developed by Zhao (1996). The highest resolution 
haloes to date are a series of eight clusters, several with 
~ 20-30 x 10 6 particles (Fukushige, Kuwai, & Makino 2004), 
which have central slopes steeper than NFW and shallower 
than M99. 

The level of "universality" of the density profile is a 
matter of debate (e.g. Tatitsiomi et al. 2004). Jing & Suto 
(2000, 2002), found central density cusps of r _1 '\ r -1 ' 3 , and 
r~ for a simulated halo with cluster, group, and galaxy 
mass, respectively. A similar range of inner slope values was 
found in a recent set of high resolution haloes (Hayashi et 
al. 2003; Navarro et al. 2004). The central cusp is especially 
sensitive to flattening caused by poor resolution or other nu- 
merical effects (e.g. Moore et al. 1998), so slightly different 
numerical techniques might produce significantly different 
density profiles, making it sometimes difficult to compare 
results of different authors. The appearance of near univer- 
sality in many previous studies could be due to the fact 
that most simulations have modelled objects in the range 
of galaxies to clusters; in this mass-range the effective slope 
of the linear power spectrum (n, where P(k) oc k n ) is n ~ 
-2, implying cusps of roughly NFW slope (Syer & White 
1998; Subramanian, Cen, & Ostriker 2000). Ricotti (2003), 
using a large set of haloes with 10 4 to 10 5 particles, found 
considerably flatter cusps for high-redshift low-mass haloes, 
with r~ ' 5 at z > 10 for 10 8 /i _1 M©. Flat, low mass haloes 
match model predictions (Syer & White 1998; Subramanian, 
Cen, & Ostriker; 2000), wherein the central slope varies as 
r (9+3n)/(5+n) . See however, Moore et al. (2001), who find a 
r" 1 ' 3 cusp in a 10 8 /i _1 M© halo at redshift four. In CDM 
models, any dependence on n would be manifested as a de- 
pendence on halo mass. In the ACDM model, n asymptoti- 
cally approaches -3 for low mass haloes. As n nears -3, M*, 
the characteristic mass of collapsing haloes as a function of 
scale factor (see section 4.3) diverges, so haloes of all masses 
collapse nearly simultaneously. If one models halo formation 
as the assembly of spherically symmetric shells of material 
whose density is largely determined by the scale factor of the 
universe, then the density profile should be shallower when 
n nears -3. 

The characteristic radius, r s , indicates the size of the 
central density region and is usually defined in terms of 
the concentration parameter, c = r v i r /r s . NFW, as well 
as a number of other authors found that the concentration 
radius decreases with halo mass, even for power law cos- 
mologies where n is constant. This agrees with other sim- 
ulations (e.g. Bullock et al. 2001a; Eke, Navarro & Stein- 
metz 2001), and also with predictions based on the halo 
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model wherein nonlinear properties are predicted starting 
from Press & Schechter (1974) theory (e.g. Huffenberger & 
Seljak 2003). The concentration dependence on mass be- 
comes weaker when n is closer to -3 (Eke, Navarro, & Stein- 
mentz 2001). 

In spherically symmetric infall models, concentration 
should increase with lower mass because lower mass haloes 
form early, when the universe is more dense, so their central 
regions are assembled with comparatively higher densities 
than their outer regions. This should give them higher char- 
acteristic central densities, or equivalently, higher concentra- 
tion parameters (e.g. Eke, Navarro, & Steinmetz 2001, and 
references therein). Eke, Navarro, & Steinmetz (2001) and 
Bullock et al. (2001a) utilise spherical infall models to pre- 
dict that halo concentration should decrease with redshift for 
a given mass for similar reasons. The spherical infall model 
also gives an intuitive prediction for trends in central slope 
at a given radius. If halo mass is increasing rapidly with 
scale factor, then a shallow central slope would be expected 
because accreting matter will all have a relativity similar 
physical density. Similar, but perhaps more physically moti- 
vated, are merger models (Syer & White 1998; Subramanian, 
Cen & Ostriker 2000), which describe the evolution of the 
halo density profile as a result of hierarchical halo formation. 
In the merger models, when smaller satellites merge with the 
main halo, their orbits will decay via dynamical friction, and 
tidal stripping will add their mass to the parent halo. The 
inner slope is then set by the dependence of the concentra- 
tion (or equivalently, the characteristic density) on satellite 
mass, and so can be predicted from n. In general, if small 
haloes formed early, as is the case if n ~ 0, they should 
be dense, which should yield steeper central slopes because 
their orbits can decay to smaller radii before being tidally 
disrupted (Syer & White 1998). 

First works on halo mass profiles suffered from poor res- 
olution that made it difficult to evaluate the central slopes 
in haloes with high concentrations (i.e. at galaxy scales). 
Also, samples with higher resolution but with limited statis- 
tics could not study trends over a large mass range, differ- 
ent cosmologies, and were unable to evaluate the intrinsic 
amount of cosmic scatter. Finally, a large dynamical range 
and high mass resolution at high redshift are necessary to 
follow the early assembly of the central part of dark matter 
haloes. This is particularly difficult for small-mass haloes 
that have higher concentrations, form earlier, and then need 
to be evolved by several internal dynamical times. A careful 
choice of the softening, force errors, and number of timesteps 
is then necessary to avoid introducing spurious numerical 
trends. 

In this work, we improve upon previous studies by mod- 
elling a large set of high resolution haloes and following their 
evolution over a wide range in mass and redshift parameter 
space. We present a set of high resolution simulations cov- 
ering 3 orders of magnitude in mass, from 2xl0 11 /i" 1 M o 
to 2x 1O 14 /i _1 M0 . This allows us to search for potential 
mass dependent trends in halo profiles. Our highest reso- 
lution haloes have 4 million and 7 million particles within 
the virial radius, respectively. Ten of our haloes are from 
CUBEHI, a single simulation of a cosmological volume at 
uniform resolution. This allows us to analyse cosmological 
scatter in profile shapes with a uniform method not subject 
to systematic uncertainties associated with differing numer- 



ical parameters. Furthermore our haloes are resolved with 
several hundreds of thousands of particles to redshift of 2 
and, in a few cases, to redshift 4 or higher. 



2 NUMERICAL TECHNIQUES 
2.1 The simulations 

We use the parallel KD (balanced binary) Tree (Bentley 
1975) gravity solver PKDGRAV (Stadel 2001; see also Wad- 
sley, Stadel, & Quinn 2004) for all of our numerical simula- 
tions. Initial conditions for the simulations are set by map- 
ping particles on to a random realisation of the mass power 
spectrum, which is extrapolated to a sufficiently high red- 
shift, z s tart, that particle overdensities are safely in the lin- 
ear regime. For higher resolution within a single halo, we use 
a "renormalized volume" technique of nested resolution re- 
gions, which has been successful in a number of cosmological 
simulations (e.g. Katz & White 1993; Ghigna et al. 1998). 
First, a low resolution cosmological simulation is completed. 
Next, a halo of interest is identified. To minimise sampling 
bias, volume-renormalized haloes are selected by mass with 
the only additional constraint that they not lie within close 
proximity (2 — 3r v i r ) to a halo of similar or larger mass. 
Then, the initial conditions routine is run again to add small 
scale scale power to a region made up of high resolution 
particles that end up within approximately two virial radii 
of the halo centre, while preserving the original large scale 
random waves. This process is iterated in mass resolution 
increments of a factor of eight until the desired resolution 
is achieved. We have verified that the high-resolution haloes 
are free from significant contamination by massive particles. 

All of the simulations model a ACDM cosmology with 
fi m =0.3 and A =0.7. We normalize the density power spec- 
trum of our initial conditions such that <jg extrapolated 
to redshift of zero is 1.0, consistent with both the cluster 
abundance (see e.g. Eke, Cole, & Frenk 1996 and references 
therein) and the WMAP normalization (e.g. , Bennett et 
al. 2003; Spergel et al. 2003). We use a Hubble constant of 
h =0.7, in units of 100 km s^ 1 Mpc -1 , and assume no tilt 
(i.e. a primordial spectral index of 1). To set the initial con- 
ditions, we use the Bardeen et al. (1986) transfer function 
with r = n m (z = 0)h. See Reed et al. (2003) for further 
details on the uniform resolution CUBEHI run. 

Our high resolution simulations are listed in Table 1. For 
the volume-renormalized runs, we list the effective particle 
number of the highest resolution region rather than the ac- 
tual particle number. Softening choices are chosen based on 
empirical studies (e.g. Moore et al. 1998; Power et al. 2003). 
Force softenings are r so ft = 5/t -1 kpc for the uniform resolu- 
tion CUBEHI run, and r so ft = 1.5% times the mean inter- 
particle spacing for the volume renormalized runs. Long 
range forces are calculated by hexadecapole expansion of the 
potentials of distant tree nodes (or "cells") that subtend an 
angle less than the cell opening angle, O, chosen to be con- 
sistent with tests by Stadel (2001). O is set to be smaller at 
high redshift, when force errors due to long range gravita- 
tional forces can have a larger contribution to total forces 
since the density field is more uniform. Ti mesteps for the 
CUBEHI run are constrained to At < 0.2-^/ r so ft/a, where 
a is the magnitude of the acceleration of a given particle, 
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Table 1. Summary of our halo sample at redshift zero. For volume renormalized runs, the mass and particle number of the central halo 
is listed. N p e g is the effective particle number based on the high resolution region for renormalized runs 



M H alo N Pjhalo N PjCff rsoftC/i^kpc) 0(z>2) 0(z<2) z start L box (ft^Mpc) 



CUBEHI 


0.7-2. lxlO 14 


0.6-1. 6xl0 6 


432 3 


5 


0.7 


0.8 


69 


50 


10 clusters 


GRP1 


4xl0 13 


7.2xl0 6 


1728 3 


0.625 


0.5 


0.7 


119 


70 


Fornax mass 


CL1 


2.1xl0 14 


4.6xl0 6 


864 3 


1.25 


0.5 


0.7 


119 


70 


Cluster 


GAL1 


2xl0 12 


0.88X10 6 


2304 3 


0.469 


0.5 


0.7 


119 


70 


Milky Way 


GRP2 


1.69X10 13 


0.38X10 6 


864 3 


1.25 


0.5 


0.7 


119 


70 


Group 


DWF1 


1. 88x10" 


0.64X10 6 


4608 3 


0.234 


0.5 


0.7 


119 


70 


2 Dwarfs 




1.93X10 11 


0.66X10 6 
















n = 


1.9X10 14 


0.54X10 6 


432 3 


2.5 


0.5 


0.7 


799 


70 


P oc k° 


n = -1 


2xl0 14 


0.55x10 s 


432 3 


2.5 


0.5 


0.7 


269 


70 


P oc k- 1 


n = -2 


1.6xl0 14 


0.45 xlO 6 


432 3 


2.5 


0.5 


0.7 


99 


70 


P oc k- 2 


n = -2.7 


2.9xl0 13 


0.82 xlO 5 


432 3 


2.5 


0.5 


0.7 


79 


70 


Pock' 27 




10' 



10 14 



= 10 12 =- 



10 11 - 



10 10 



Number particles 
Within Rvir 
O > 4xl0 6 
■ 2-4X10 6 
l-2xl0 6 
A 0.5-lxlO 6 
» 0.1-0.5x10" 



3 3 
Redshift 



Figure 1. Our set of haloes for density profile analyses, including 
lower resolution tests. 



for the CUBEHI run, and At < O.l75^r so[t /a for all other 
runs. These timesteps are consistent with convergence tests 
for variable timestep runs by Power et al. (2003), where 
0.2-^/ r so f t I a is found to be sufficient for the central 

is 1 



At 



regions of haloes. The number of periodic replicas, n r 
for all simulations; n r determines the number of copies of 
the box for gravity calculations, and hence the accuracy of 
the periodic force. Starting redshift, r so ft, O, and n r are all 
tested for the CUBEHI run in Reed et al. (2003). 

For our volume-renormalized galaxy, group, and cluster 
runs, we have analysed lower resolution runs which we use 
for a numerical resolution study. In the interest of limiting 
our study to only the highest resolution haloes, we only cal- 
culate density profiles for the ten most massive haloes in our 
CUBEHI run, each of which has ~10 6 particles. For similar 
reasons we only follow the evolution of haloes to redshifts 
where their particle number exceeds 10 . In Fig. Q we plot 
the properties of the full sample of haloes that we analyse, 
including our low resolution versions of the renormalized 



runs GAL1, GRP1, and CL1 when they exceed 10 5 parti- 
cles. The n — 0, —1, —2, & — 2.7 runs are used to study the 
effects of power index spectral slope on halo structure. The 
initial power spectrum for these runs is given by P oc k n 
normalized to the same as as all other runs. These power 
spectra are based on the same random waves as run CL1. 



2.2 The analysis 

Our virialized haloes are selected with the (SO) algorithm 
(Lacey & Cole 1994) utilising the spherical tophat model 
of Eke, Cole, & Frenk (1996) in which the ACDM virial 
overdensity, A v i r , in units of critical density is approximately 
100. To follow the evolution of an individual halo, we "mark" 
a few hundred particles at the density peak of the halo at 
z = and trace those particles back to higher redshifts, 
when they are in the core of the largest progenitor haloes. 
When there are multiple progenitors, we use the progenitor 
with the deepest potential. 

To calculate density profiles without excessive particle 
noise, we developed a novel kernel-based algorithm 1 (Merritt 
& Tremblay 1994); see Appendix for a full description. Both 
the width and shape of the kernel are varied with radius; 
the variation in shape is significant near the origin, where 
a symmetric kernel would "overflow" the r = boundary. 
The window width must be carefully chosen to reduce Pois- 
son noise ("variance") without oversmoothing ("biasing") 
the profile. In general, it may be shown (e.g. Scott 1992, 
p. 130) that when the window width is chosen to minimise 
the mean square error of the estimate, most of the error 
will come from the variance. Window widths large enough 
to eliminate the "wiggles" will generally bias the slope. In 
addition, the window width should vary with local particle 
density (Abramson 1982), roughly as p" 1 ^ 2 . We used a ker- 
nel window width that varied as r ' 5 set at ho.i = .005r V ir at 
0.1rvi r as it yields profiles and profile slopes in good agree- 
ment with binned profiles created with TIPSY 2 ), and it pre- 
serves the central cusp and major substructure. 



1 Code available at |htt p : / / www .rit ,edu/~ drmsps/ inverse . 

2 TIPSY is available from the University of Washington N-body 
group: http:/ /hpcc. astro. washington.edu 
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Figure 2. Top panels: Density profiles of our galaxy halo (GAL1, left) and our highest resolution halo (GRP1, right) at multiple 
resolutions. Long arrows show resolution criteria from Moore et al. (1998); medium arrows show Power et al. (2003) criteria; and short 
arrows show Fukushige & Makino (2001) criteria. Softening is indicated by vertical lines at top of plot windows. Profiles are constructed 
via a kernel estimator (Appendix). Bottom panels: Derivative of the density profile of above haloes at the same resolutions. Density 
slopes are calculated using a rolling average of Ar/r v ; r ~ ±30%. 



3 RESOLUTION CRITERIA 

In this section, we utilise lower resolution versions of our 
renormalized- volume simulations to examine several resolu- 
tion criteria. Previous authors have proposed empirical and 
theoretical criteria to define the minimum radius at which 
the density profile can be considered to be "resolved" (see 
Diemand et al. 2004 and references therein, whose discus- 
sion we summarise here). The main numerical issue is the 
discreteness caused by the fact that N-body particles are ex- 
tremely massive compared to dark matter candidates in the 
"real" universe. This discreteness means that particles will 
undergo two body interactions, which change their velocity 
by a significant amount. Two body relaxation effects van- 
ish as N p approaches infinity. The two body relation time, 
defined as the time it takes for particle energy to change 



by order unity, is shortest near halo centres where density 
is high. In CDM haloes, two body interactions tend to add 
energy to the low velocity cusp particles, so two body re- 
laxation has the effect of flattening the inner cores of sim- 
ulated haloes. After several Hubble times, haloes become 
nearly isothermal, and the energy transport reverses direc- 
tion. Power et al. (2003) and Fukushige & Makino (2001) 
have considered the relaxation rate at z — and found that 
haloes are resolved down to radii where the relaxation time 
is equal to the Hubble time and three Hubble times, respec- 
tively. Moore et al. (1998) and Ghigna et al. (2000) offer 
an empirical fit finding that the minimum resolved radius 
is r min ~ / m ean/2 where l mea n = (47r/3) 1/3 A r ,7 1/ ' 3 , based 
on simulations of identical haloes at different resolutions; 
see Splinter et al. (1998) and references therein for discus- 
sion relating resolution to mean particle spacing. Because 
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most two body relaxation occurs early, when particles are 
in small haloes, particles in higher-mass haloes can suffer 
significantly more relaxation (Diemand et al. 2004). This is 
due to the later formation time of massive haloes, which 
means that their particles have spent more time in small jV p 
progenitors where two-body relaxation is larger (Diemand 
et al. 2004). We therefore test our resolution criteria over a 
range of masses. 

We use identical haloes of varying resolution to empiri- 
cally identify radii we consider to be well resolved and com- 
pare our results with other resolution studies. Fig. [2] shows 
the density profiles of our galaxy and group for which we 
have multiple resolutions and indicates the resolution crite- 
ria from the above studies. We have made similar resolu- 
tion studies of cluster CL1. In each case, none of the three 
resolution under-estimates the radii of divergence between 
the point of divergence between lower and higher resolu- 
tion haloes over the full mass range. However, the Power et 
al. and Fukushige & Makino criteria are more conservative 
for our lower resolution haloes, as they scale more steeply 
with particle number. Fig. [5] also shows the slopes of the 
density profiles, dlnp/dln(r/r v i r ), for the same haloes. To 
reduce noise, profile slope is calculated based on a rolling 
average of the kernel-based density profile of roughly ±30% 
in radius 3 This slope estimation is sufficient for our pur- 
poses; however, we note that an optimal method of obtain- 
ing low-noise profile slopes is to compute the density deriva- 
tive directly from a kernel based density profile that was 
made using a larger kernel window. The profile slopes seem 
more sensitive to particle resolution, and thus appear to be 
accurate down to minimum resolved radii that are ~ 50% 
larger than inferred from the density profiles. At low particle 
numbers, both Power et al. and sometimes the Fukushige 
& Makino criteria appear to be over conservative, which 
suggests that because of their steeper particle number de- 
pendence, these criteria may not be conservative enough for 
haloes with very large N p . All of our haloes seem well re- 
solved down to a radius a little larger than i mea n/2- We thus 
utilise r m i n = Np 1//3 , which is 25% larger than Z mcan /2- This 
empirical criterion seems to best match the dependence of 
r m in on particle number in our simulations. We note that 
one should not expect this criteria to be valid for haloes 
with vastly different central densities or particle numbers 
than modelled here, as the relation between resolved ra- 
dius should depend on the central halo density and other 
physical properties in addition to particle number. We have 
performed similar resolution tests of z = 1 outputs, where it 
appears that most haloes are resolved to slightly better than 

— 1/3 

N p , probably because particles have had less time to un- 
dergo two body interactions. At this high redshift, the Power 
et al. formula still gives a conservative resolution limit, but 
the Fukushige & Makino criteria becomes less conservative 
than Z mean /2 for our highest particle numbers. In sum, the 



3 Since the kernel density estimate is a continuous function of ra- 
dius, a better way to compute the derivatives would have been via 
analytical differentiation of i>(r), using a somewhat larger window 
width to compensate for the increased variance of the derivative 
as compared with the function itself (e.g. Scott 1992, p. 131). We 
recommend that this procedure be followed in the future, though 
it has no significant effects on our results. 



minimum resolved radius of the density profile is well de- 
scribed by N p for a wide range of halo masses and red- 
shifts, so we adopt this criteria for the remainder of the 
paper. 



4 RESULTS 

4.1 The full halo sample 

In Fig. [3] we present a plot of the density profile of all of our 
haloes at redshifts of zero and one. Large substructure is ap- 
parent in the profiles which have been constructed with the 
kernel algorithm described earlier. Each halo has a unique 
density profile, even in the inner regions where there is lit- 
tle substructure. The cusp size and slope vary from halo to 
halo. The redshift one profiles are normalized in terms of 
the redshift zero virial radius and critical density in physi- 
cal (non-comoving) units, and are plotted out to the z = 1 
virial radius. There is little evolution of the profile between 
redshift zero and one. In Fig. 0] we plot the z = slopes of 
the same density profiles, dlnp/dln(r/r v j r ), which is calcu- 
lated based on a rolling fit of approximately ±30% in radius. 
Substructure is prominent in the outer regions of each halo. 
Only a few of the halo density profiles appear to converge to 
an asymptotic inner slope; the rest continually flatten all the 
way down to the innermost resolved radii, though at a rate 
generally consistent with an asymptotic slope parameter of 
r _1 or steeper. The NFW and M99 curves are plotted for 
the best fit concentration parameters. The innermost slope 
at the inner resolution limit ranges between r -1 ' 1 and r" 1 ' 7 , 
which implies that the halo to halo cosmic scatter in slope 
is approximately bounded by the NFW and M99 profiles. 

In Fig. the z — halo density slopes are plotted 
versus a two parameter fit where the inner asymptotic slope 
parameter is allowed to vary in addition to the concentration 
parameter. The density is thus given by: 



' (c 7 r/r vir )^[l±(c 7 r/r vir )] 3 -^' 

where 7 is the asymptotic inner slope parameter and c 7 is 
the concentration parameter obtained when 7 is a free pa- 
rameter. This two parameter fit produces a visually better 
fit to the density profile in most cases. The range in (-)7 for 
the 2 = haloes is -1 to -1.7. Note that there is a partial de- 
generacy between 7 and the concentration radius. A pseudo 
(because the data points are correlated) x 2 per degree of 
freedom for each of these fits shows substantial improvement 
over NFW or M99 fits; see Table 2. The profile fits are based 
on a least squares method, where a Poisson uncertainty is 
estimated for each point based on the effective number of 
particles for the density given by the kernel-based profiles, 
with logarithmically spaced bins. 

4.2 Redshift Evolution 

In order to gain an understanding of the physical effects 
that set the inner slope of halo density profiles, we plot the 
evolution of the profile slope in terms of the z — virial ra- 
dius in non-comoving coordinates for our four best haloes. 
This allows us to ignore the effects that expansion of the 
universe or evolving r v i r have on the power law slope of 
the profile. We see in Fig. |S| that the density profile inner 
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Table 2. Our best-fit concentration parameters for the NFW profile (column 3) and the M99 profile (column 
5) followed by their respective pseudo-x 2 per degree of freedom goodness of fits at redshift 0. Columns 7 and 
8 are for a two parameter fit where 7 is the central slope parameter, followed by its pseudo-x 2 m column 9. 
Column 10 is the characteristic radius for this two parameter fit 
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Figure 3. Density profiles of our 16 haloes sample for at z = Figure 4. Density profile slopes of our 16 haloes sample for at 

(solid) and z = 1 (dashed). Here we rescale the 2 = 1 profiles 2 = 0. The best- fit concentration NFW (red, long dashed) and 

in physical (non-comoving) coordinates normalized to the critical M99 (blue, short dashed) profile slopes are plotted. Profiles are 

density and virial radii at z = such that a value of r/r virj0 corre- plotted to the minimum resolution criteria of N~ 1/3 . Slope is 

sponds to the same non-comoving distance from the halo centre calculated based on a rolling fit of ±30% radial width, 
at all redshifts. Profiles are plotted to the minimum resolution 
criteria of 7V p . 

profile slope shows no significant change. The large fluctua- 
tion in the profile slope at z — 3 for the galaxy halo and at 

slope for each halo evolves very little in non-comoving co- z = 5 for the dwarf halo are most likely due to the presence 

ordinates, plotted as dlnp/dln.(r/rvi r ,o), a result also found of subhaloes that are disrupted earlier. Whatever physical 

by Fukushige & Makino (2001) and Fukushige et al. (2004). mechanism is responsible for setting the density profile, it 

In our cluster CL1, the slope slowly steepens with time at must have largely occurred at very high redshift. 

z 2, but in our other less massive haloes, the inner density The lack of evolution in the physical densities means 
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0.01 0.1 0.01 0.1 



r / r vir,0 

Figure 6. Evolution of the slope of the density profiles of four haloes plotted against r v i r at z = in non-comoving (physical) coordinates 
as in Fig. |3 Mass beyond the virial radius of each epoch is ignored. Note that if plotted simply in terms of r/r v ir the haloes would 
appear to flatten with increased redshift. 



that the profile of an individual halo is substantially shal- 
lower at high redshift in terms of r/r v i r . Such apparent steep- 
ening of the slope with time is merely a scaling issue due to 
the growth of the virial radius as mass is added to the outer 
regions of the halo. Here we note that the profile slope is 
almost never shallower than the NFW asymptotic value of 
r _1 . Only in the very inner region of our largest cluster CL1 
at 2 = 3, where the slope trends toward r~ 0,75 at the in- 
nermost resolved radius (not plotted), is there a slight hint 
that NFW slope may be too steep at high redshifts, but this 
could be simply a result of a subhalo just beyond that ra- 
dius creating a local density minima, or it could be due to 
artificial numerical resolution effects at high redshift where 
the halo resolution is lower. 



4.3 Trends in Profile Concentration and Inner 
Slope 

The profile concentration and thus the inner slope at a 
given radius are predicted by NFW and others to be a 
function of mass and redshift, when considered in terms 
of r/rvir, with characteristic radius r s increasing with in- 
creasing M/M*. Here, M* is the characteristic mass of col- 
lapsing haloes defined by the scale at which the rms linear 
density fluctuation equals the threshold for non-linear col- 
lapse (i.e. a(M*(z)) = S c ). We have measured the NFW 
concentration parameter, cnpw, by forcing our profiles to 
an NFW profile, and performing a least squares fit. In Fig. 
[7| we plot the concentration parameter for our set of haloes, 
and we show the Eke, Navarro, & Steinmetz (2001) predic- 
tion for cnfw- The concentration dependence on M/M* for 
our haloes is significantly steeper than predicted by NFW for 
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Figure 5. Density profile slopes of our 16 haloes sample at z = 0. 
The long dashed curve (red) shows the results of a two parameter 
density profile fit in which both the inner slope parameter, 7, and 
the concentration parameter, c 7 , are allowed to vary. The top 
number in each window is the best-fit c 7 , and the bottom is the 
value of the best-fit (-)7- 



M < M» , and within the scatter of our haloes for M > M* . 
Measured values of the inner slopes also show clear trends 
with mass and redshift, shown in Fig. |H] though this is due at 
least in part, to the degeneracy between slope and concentra- 
tion. Here, we have measured the average value of the power 
law slope dlnp/dln(r/r v i r ) between 2% and 5%r v i r . The 
asymptotic inner slope parameter, 7, from a two-parameter 
fit of 7 and concentration as given by equation 3, is plot- 
ted in Fig.|U] and has a weak dependence on M/M* for our 
haloes, with a large scatter. The median value of 7 in our 
sample trends toward shallower inner slopes with increasing 
M/M», given by: 

7~1.4-0.08Log 10 (M/M.), (4) 

with a scatter of A7 ~ ±0.3 for our haloes, and valid for 
haloes of 0.01M, to 1000M,. Fig. [T01 shows that the con- 
centration parameter in this two-parameter fit, c 7 , shows a 
significant trend toward higher values as M/M, decreases. 
The trend is weaker and has significantly more scatter than 
the forced NFW fit cnfw dependence on M/M* (see Fig. 
Q . Adopting a power-law parameterisation as in e.g. Huf- 
fenberger & Seljak (2003), the median c 7 for our haloes is: 

c 7 ~8.(M/M,)- al5 , (5) 

with a M/M, dependent scatter roughly equal to ±c 7 . 
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Figure 7. Best fit NFW concentration parameter r v i r /r s for our Figure 9. The asymptotic inner slope parameter, 7, from a two 

set of haloes as a function of M/M*. Empirical predictions by parameter fit of our 16 halo sample. 7 and the concentration are 

Eke, Navarro, & Steinmetz (2001) are given by open circles. both allowed to vary. Solid line fit is given by Eg. El 
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Figure 8. Density profile slopes of our 16 haloes sample averaged 
over the range of 2% to 5% r v i r . Haloes not resolved to 2% r v i r 
are not plotted. 



Figure 10. The concentration parameter from a two parameter 
fit of our 16 halo sample, c 7 . The asymptotic inner slope param- 
eter, 7, and the concentration are both allowed to vary. Solid line 
fit is given by Eg. IS1 
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4.4 Cosmological Variance and Stability of 
Profiles 

In our CUBEHI simulation, we examine the stability in the 
profiles of our set of ten clusters over timescales separated 
by Az ~ 0.01 for z ^ 0.2; Fig. 1111 In the outer regions, 
orbiting substructure creates substantial scatter, however, 
the inner density profile is relatively stable. Thus, differences 
between the central density profiles of haloes of similar mass 
reflect different inherent properties of the halo, rather than 
temporal effects of orbiting substructure. 

In the hopes of understanding what physical processes 
are responsible for setting the central density profile of each 
halo, we examine the evolution of the main halo and its 
progenitors. We construct merger histories for each cluster 
and examine a number of properties, including: evolution of 
the central cusp mass concentration; cluster accretion his- 
tory; total collapsed progenitor mass; and angular momen- 
tum profiles. To follow the mass accretion history of the 
cluster and its progenitors, we use the friends-of-friends al- 
gorithm (FOF, Davis et al. 1985). In order to follow the evo- 
lution of the mass concentration that makes up the cluster 
cusp, we use SKID 4 (Stadel 2001), which is able to identify 
bound mass concentrations independently of environment. 

In Fig. 1121 we plot the mass evolution of the central 
SKID progenitor. The two haloes with the earliest forming 
central SKID halo are also those that have density profiles 
most closely matching the M99 profile (cube4 and cube6; see 
Fig- El Table 2), both with slopes that steepen slowly toward 
r -1,5 , beginning at radii of roughly 10%r\,i r , implying large 
concentration radii. If not simply a coincidence, then this 
implies that the central cusp material is assembled earlier in 
haloes with steeper central slope parameters. We then con- 
sider the effects of the mass accretion history of the cluster 
on the final density profile, and find that accretion history 
correlates with halo concentration, not with cusp slope. Fig. 
1131 shows that the three haloes with the highest concentra- 
tion undergo a phase of rapid growth at z ~ 2 — 8, making 
their normalised mass temporarily ~ 3 times larger than 
the other seven clusters, which experience nearly uniform 
accretion rates. We have also examined the evolution of the 
total collapsed progenitor mass, and find a similar correla- 
tion with cluster concentration. At redshifts of ~ 10, the 
same three highly concentrated haloes have ~ 3 times more 
total mass in collapsed progenitors (not plotted), where we 
have only considered progenitors of mass greater than 0.01 
percent of the final cluster mass. The correlation of the con- 
centration parameter with halo and progenitor collapse time 
is not surprising; in fact, as have shown in section 4.3, lower 
mass haloes, which are assembled earlier in a hierarchical 
model, have smaller concentration radii. The concentration 
trends also qualitatively agree with correlations of forma- 
tion epoch and concentration found in numerical simula- 
tions by Wechsler et al. (2002). However, a cautionary note 
is needed here. Even though in the CUBEHI simulation, the 
halo masses differ by only a factor of ~ 3 and thus cover 
a narrow range in median concentration parameter (see eq. 
5), the two least massive and hence most poorly resolved 
clusters, are also the two most concentrated. A larger set 
of higher resolution haloes is needed to conclusively rule out 

4 SKID available at: |http~/ hpcc. astro. Washingt on. edu| 
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the possibility that haloes resolved with fewer particles (and 
identical softening) lead to higher concentrations. In Fig. 1141 
we plot the angular momentum parameter A = j/\/2v c r, 
from Bullock et al. (2001b), where j is the specific angular 
momentum. We see no correlation of profile concentration or 
central slope with A in either the inner or the outer regions. 
We also have examined the angular momentum profiles at 
high redshift, but find no clear correlation with halo concen- 
tration or central slope. This is puzzling given that merger 
histories, which correlate with the density profiles, should 
also correlate with the angular momentum distribution. 

4.5 Power Law Cosmologies 

In order to understand the effects of the power spectral slope 
index, n, we have simulated the renormalized volume cluster 
CL1 in scale free cosmologies with a range of values for n. 
Here the initial density fluctuation power spectrum is given 
by P oc k n , and is normalized to erg = 1.0 with S7 m = 0.3 and 
A = 0.7. Fig. 1151 shows the density profile for a cluster with 
power law initial conditions given by n =0, -1, -2, and -2.7, 
followed by the corresponding density profile slopes in Fig. 
1161 Here we also include plot the Power et al. density criteria 
as we have not done convergence tests specifically for P oc k n 
cosmologies. Note that the n = —2.7 run has significantly 
less mass because of the ag = 1.0 normalization. There is a 
clear trend that steeper power spectra yield density profiles 
with flatter slopes and larger concentration radii. Visually, 
the n = cluster displays much more prominent substruc- 
ture, due to its proportionally stronger small scale power, 
with many nearly spherical haloes and few filaments. A two 
parameter profile fit yields 7 = 1.5 for n — and 7 = 1.1 
for n — —2.7 (see Table 2). A number of authors (Hoffman 
& Shaham 1985; Crone, Evrard, & Richstone 1994; Cole & 
Lacey 1996; Syer & White 1998; Subramanian et al. 2000; 
Huffenberger & Seljak 2003; NFW 1997; Eke, Navarro & 
Steinmetz 2001; Ricotti 2002) have suggested a power law 
dependence of the density profile that qualitatively agrees 
with our power law simulations. However, Syer and White 
(1998) predict that an n — —2.7 power law cosmology should 
have an inner slope of r -0,4 , which is much shallower than 
seen in our haloes, though not necessarily inconsistent if the 
slope flattens only at very small radii. The same model pre- 
dicts an inner slope of r~ 8 for the n — cosmology, which 
is very close to the inner slope at the minimum resolved 
radius of our n — halo. 



5 DISCUSSION AND CONCLUSIONS 

(1) Density profile trends and the CDM "cusp problem". We 
follow the evolution of sixteen haloes over a large range 
in parameter space, with masses corresponding to dwarfs 
through clusters back to z = 3 for our highest resolution 
cluster and 2 = 5 for two dwarfs. By using identical haloes 
at varying resolution, we have shown that our haloes are 
likely to be free from biases related to numerical resolu- 
tion. None of our simulated haloes has a density profile slope 
significantly shallower than r _1 down to the minimum re- 
solved radii of 0.5% to 2%. These steep cusps are similar 
to those found in previous CDM simulations, and appear 
to be in conflict with reported rotation curves from dwarfs 
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Figure 11. Density profile slopes of our 10 clusters from the CUBEHI simulation. Here we plot the density profile at intervals separated 
by Az~0.01 up to z =0.2. Note that these profiles are binned, rather than kernel-based. From top to bottom and left to right, the haloes 
are plotted based on descending mass. The lower right plot is the average of the slope value for all ten clusters. The NFW profile with 
the best fit concentration value is denoted by the solid (red) curve for each halo 



and LSBs. Our resolution experiments confirm that steep 
cusps are formed regardless of numerical resolution, down 
to the minimum resolved radius of each halo. Previous high 
resolution simulations discussed earlier have mostly mod- 
elled haloes with masses of lO 12 fr -1 M0 or higher, leaving 
the possibility that a mass dependence on the density pro- 
file might be able to solve the cusp problem. However, our 
results show that the observed conflict with ACDM density 
profiles gets worse when considering simulated dwarfs. The 
enclosed mass at a given radius near the halo core is higher 
for simulated haloes of lower mass, since the concentration 
radius is smaller and the measured slope is steeper. 

(2) Convergence of profile slopes. We confirm that the 
physical slope of haloes remains stable over time. The in- 
ner regions of haloes appear to be composed of mass as- 
sembled at very high redshift, implying that present-epoch 
density-profiles are determined by the high redshift merger 
history. We do see some apparent evolution in the physical 
slope at very high redshifts for our highest resolution clus- 
ter and group. This suggests that for group and clusters-size 
haloes, events occurring at z ~ 2 — 4 are partially responsible 
for determining the final shape of the density profile, while 
galaxies and smaller mass haloes have their inner density 
profile shape almost entirely determined at higher redshifts. 



The progenitor region of the host halo, because of its large 
scale density enhancement, is the likely site of early-forming 
small haloes, many of which will merge to form the cusp. Af- 
ter the universe has expanded by a few factors since the cusp 
material is assembled, the characteristic densities of haloes 
merging into the main halo should be lower than typical 
densities within the main cusp, because of their generally 
later assembly epoch. Consequently, merging subhaloes will 
likely be disrupted before dynamical friction can bring them 
near the centre of the main halo, where they can affect the 
central density profile. In this scenario, halo density profiles 
could converge to roughly flat profiles (in terms of r/r v i r ) at 
extremely high redshifts. CDM haloes at low redshift might 
then be expected to have a very small flat inner core, but at 
radii much less than the inner ~ 1% that can be currently 
probed by simulations or observations. 

(3) Lessons from the power law cosmologies. The de- 
pendence of the concentration parameter on M/M* is most 
likely the result of a mechanism analogous to what sets the 
dependence of concentration on the spectral slope in our 
power law cosmologies. In the power law cosmologies, the 
ratio of local small-scale to large-scale power depends pri- 
marily on the spectral index. When there is lots of small 
scale power, the central density profile is generally steeper 
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Figure 12. Evolution of the mass of the central SKID progenitor halo for the same ten clusters as Fig, 1111 For reference, the most 
massive halo is replotted in green (dashed curve) in each plot window. 



at any specific radius than when there is little small scale 
power. In the case of a shallow spectral index, subhaloes 
form early, which implies that they have small characteristic 
radii and thus large characteristic densities, whereas steep 
spectral indices have late-forming subhaloes. Similar forma- 
tion epoch arguments can also explain the steeper profiles of 
small M/M, haloes. For smaller M/M, haloes, many sub- 
haloes should form earlier relative to their hosts, since the 
subhaloes lie at lower a fluctuations in the density field rel- 
ative to the local fluctuation on the scale of the host halo. 
Alternatively, smaller M/M, haloes, which generally form 
earlier, and hence when the mean density of the universe is 
evolving more rapidly, will have subhaloes and cusp mate- 
rial assembled over a wider range of densities than for large 
M/M, haloes. Haloes with small M/M, are thus similar 
to haloes with a shallow spectral index, because they have 
subhaloes that form earlier at higher densities relative to the 
host halo. In models where the inner slope is formed from 
disrupted subhaloes infalling via dynamical friction, lots of 
small-scale power should give the main halo a larger con- 
centration parameter. In fact, when there are lots of dense 
subhaloes more of them will be able to reach closer to the 
core, steepening the halo nearer to the centre. 

(4) Is there an asymptotic cusp? A strength of the NFW 
profile, is that varying just the concentration radius can yield 
reasonable fits to haloes with very different inner slopes at 



scales currently resolved by simulations. However, there are 
no compelling theoretical arguments that the halo profile 
converges to a slope of r _1 at smaller radii. The analyt- 
ical models mentioned in the introduction imply that the 
core slope should be a function of M/M, or power spec- 
tral slope. It is possible, and theoretically motivated, that 
haloes of higher M/M* or steeper spectral slope than what 
we have simulated here would actually have a shallower slope 
than r _1 , which would be inconsistent with the NFW pro- 
file. Most of our haloes have density profiles that get ever 
shallower with decreasing radius down to our minimum re- 
solved radius, although they are still consistent with asymp- 
totic cusp slopes of r _1 or steeper. Additionally, a few of our 
haloes appear to have converged to slopes at r -1,5 or steeper, 
which would be inconsistent with the NFW profile, though 
because of the degeneracy between central slope and concen- 
tration parameter, the NFW profile is not definitively ruled 
out. Simulations of haloes with very high M /M, , very steep 
power spectral indices, or orders of magnitude more parti- 
cles able to probe much farther inward, should be able to 
test whether or not the r" 1 NFW central slope corresponds 
to a minimum possible slope set by the physical processes 
of CDM halo formation. Navarro et al. (2004) recently pro- 
posed a new profile form based on fits to a set of high reso- 
lution haloes. The Navarro et al. profile differs from ours in 
that it has no asymptotic cusp, but instead continually fiat- 
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tens with decreasing radius. It does not become shallower 
than r _1 until radii smaller than currently resolvable, so it 
is generally consistent with most of our haloes. However, our 
profile form is more flexible in that it is able to better match 
haloes that have steep asymptotic cusps with large concen- 
tration radii, and is still a good match to those haloes that 
continuously flatten down to the minimum resolved radii 
in our simulations. Future simulations will likely need to 
probe below ~0.00I r v i r in order to determine whether dark 
haloes have an asymptotic central cusp. The mass distribu- 
tion at such small radii has a strong effect on the flux of 
hypothesized dark matter annihilation signals that may be 
detectable via 7-rays from the galactic centre (e.g. Stoehr et 
al. 2003; Evans et al. 2004). 

(5) Profile scatter and the two parameter profile When 
we plot our haloes with the best fit concentration param- 
eters, neither the NFW nor the M99 function provide a 
good fit to all of the haloes. However, because each halo 
has a unique profile determined by poorly understood and 
complex high redshift events, any possible single parameter 
profile will be unlikely to fully describe all haloes. Our two 
parameter fit describes the general trends with M/M* , pro- 
viding a significant improvement over NFW and M99, but 
it still does not account for the large halo to halo scatter 
at a given mass and redshift. Our evidence in the CUBEHI 
clusters for a likely correlation between the cusp material 



assembly epoch and cusp slope, and between halo accretion 
history and concentration radius, suggests that the density 
profiles might be at least partly determined from the evo- 
lutionary history of the progenitor haloes and substructures 
clumps. Larger sets of high resolution haloes should be able 
to quantify the role of mergers and other stochastic processes 
in shaping the density profile. 



ACKNOWLEDGEMENTS 

We are grateful to the anonymous referee for helpful and 
constructive comments and suggestions. We thank Lucio 
Mayer for assistance with one of the runs. DR has been 
supported by the NASA Graduate Student Researchers Pro- 
gram. DR was partially supported by PPARC. FG is a David 
E. Brooks Research Fellow. FG was partially supported by 
NSF grant AST-0098557 at the University of Washington. 
TRQ was partially supported by the National Science Foun- 
dation. Simulations were performed on the Origin 2000 at 
NCSA and NASA Ames, the IBM SP4 at the Arctic Region 
Supercomputing Center (ARSC) and at CINECA (Bologna, 
Italy), the NASA Goddard HP/Compaq SC 45, and at 
the Pittsburgh Supercomputing Center. We thank Chance 
Reschke for dedicated support of our computing resources, 
much of which were graciously donated by Intel. 



© 2003 RAS, MNRAS 000.HHT71 



Density Profiles 15 




Figure 14. Angular momentum profile of the ten CUBEHI clusters shown in the previous three figures. 
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APPENDIX A: KERNEL ROUTINE FOR 
DENSITY PROFILE ESTIMATION 

Here we present the algorithms which we used to derive 
smooth estimates, 0(r) and E(-R), of the particle number 
density and surface density profiles from the iV-body posi- 
tions. 

The routines in MAPEL (Merritt 1994) allow one to de- 
rive maximally unbiased estimates of v and E using penalty 
functions that embody the approximate power-law nature of 
these functions. However the MAPEL routines are relatively 
slow, and this fact presented difficulties when constructing 
estimates using the N ~ 10 6 particle data sets. Kernel based 
algorithms are faster but potentially more biased; however 



http 
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we found that density profiles produced with MAPEL and 
kernel methods are in excellent agreement down to our min- 
inum resolved radii. Thus we have used the kernel estimator 
in our analyses. 

Our derivation follows that in Merritt & Tremblay 
(1994). In the absence of any symmetries in the particle 
distribution, a valid estimate of the number density v cor- 
responding to a set of particle positions is 



\ " 



— K 



h 



(Al) 



where h is the window width and K is a normalized kernel, 
e.g. the Gaussian kernel (see Fig. Al): 

1 



K{y) 



-K72 



(A2) 



(2tt) 3 / 2 

Now imagine that each particle is smeared uniformly 
around the surface of the sphere whose radius is r; and whose 
origin is at (0,0,0). If the density profile is actually spher- 
ically symmetric, this smearing will leave the density un- 
changed; if not, it will produce a spherically symmetric ap- 
proximation to the true profile. The spherically-symmetrized 
density estimate is 



l - 3 l-j^j d e^eK{^, (as 



|r - r { | 

r 2 + r 2 — 2rTi cos I 



a) 



(A3b) 
(A3c) 



where 6 is defined (arbitrarily) from the ivaxis. This 
may be written in terms of the angle-averaged kernel K: 



K(r,n,h) 
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— I d<j> I d9 sin<9 K 



2rn cos 6 



Substituting for the Gaussian kernel, we find 

x e~ ir * +r2)/2h2 smhim/h 2 ). 
A better form for numerical computation is 



(A4a) 

(A4b) 
H ) (A4c) 



(A4) 



K{r,n,h) = 



2(2tt) 
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We want to vary the window width with position in 
such a way that the bias-to-variance ratio of the estimate 
is relatively constant. Let hi be the window width associ- 
ated with the ith particle. The density estimate based on a 
variable window width is 



JV 

i — 1 % 



(A6) 



The optimal way to vary hi is according to Abramson's 
(1982) rule: 



hi cuv a (n), 



= 1/2. 



(A7) 



Since we don't know v(r) a priori, we instead varied hi ac- 
cording to 

h.xr . (A8) 

We found that f3 = 1/2 gave good results, which is reason- 
able since density profiles are close to v ~ r . One could 
improve on this by first constructing a pilot estimate of v 
then using Abramson's rule. 

The surface density profile could be computed via sim- 
ple projection of 0{r). Instead, we computed E(-R) directly 
from the coordinates projected along one axis. The two- 
dimensional kernel estimate of E(R) in the absence of any 
symmetries is 



N 

!=1 



jj|a-R«l 



(A9) 



where K' is the two-dimensional Gaussian kernel, 

K 'ly) = J-e-y 2 / 2 . 
yy ' 2tt 

Now smear each particle uniformly in angle <j> at fixed Ri. 
The density estimate becomes 



(A10) 
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K> (|) # , 



d = Rt+R - 2RRi cos <t>. 
In terms of the angle-averaged kernel K': 
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where the last expression was derived using the Gaus- 
sian kernel; 7q is the modified Bessel function. 
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Figure Al. Kernels for the radial density estimation prob- 
lem. Window width is ft = 0.1 and particles are located at 
r= (0.1,0.2,0.3,0.5,0.7). 
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